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Abstract. We present detailed gas-phase chemical models for the envelope of the low-mass star-forming region 
IRAS 16293-2422. By considering both time- and space-dependent chemistry, these models are used to study 
' both the physical structure proposed by Schoier et al. Il2002t . as well as the chemical evolution of this region. A 

^Nj ' new feature of our study is the use of a detailed, self-consistent radiative transfer model to translate the model 

, abundances into line strengths and compare them directly with observations of a total of 76 transitions for 18 

04 ' chemical species, and their isotopes. The model can reproduce many of the line strengths observed within 50%. 

^ I The best fit is for times in the range of 3 x 10'' — 3 x 10'* yrs and requires only minor modifications to our model for 

the high-mass star-forming region AFGL 2591. The ionization rate for the source may be higher than previously 
expected - either due to an enhanced cosmic-ray ionization rate, or, more probably, to the presence of X-ray 
induced ionization from the center. A significant fraction of the CO is found to desorb in the temperature range 
of 15-40 K; below this temperature ~90% or more of the CO is frozen out. The inability of the model to explain 
' the HCS^, C2H, and OCS abundances suggests the importance of further laboratory studies of basic reaction 

, rates. Finally, predictions of the abundances and spatial distributions of other species which could be observed 

' by future facilities (e.g. Herschel-HIFI, SOFIA, millimeter arrays) are provided. 

O ■ 
^' 

^ 1. Introduction spatial distribution of the submillimeter continuum emis- 

• • . sion from dust and its spectral energy distribution (SED) 

>■ , The distribution and composition of dust and gas around gj^j^^^y ^1 j0rgensen et al. Schoier 

■ isolated low-mass young stellar objects (YSOs) is receiv- al.a002J- Through continuum radiative transfer calcu- 

^ , ing increased attention both observationally and theoret- ^^^-^^^^ ^^^^ ^^^^-^^ p^^^^l^ ^ ^ ^-p ^ ^^^^^^^^ 

CA ' ically. While the general process of low mass star forma- ri- ^ -i-i, ^ -i- -1- t^/\ 
. . T 1 ° rrrmm ^t radius r and the dust temperature structure J dust i?") 
tion is relatively well understood (see e.g. Shu et aLUHSl ean be determined self-consistently. For the gas, two ap- 
Evans[TO| and others) , many details on the chemical and p^^aches can subsequently be taken. In atmospheric chem- 
physical structure at different stages of evolution remain .^^^^^ ^^^^^ ^^^^^ commonly known as the 'forward' 
uncertain. In particular, the warm and dense gas in the 'backward' or 'retrieval' methods. In the 'empirical 
very interior of these star-forming regions provides a rich ^^^^p (..gt.ie^^p)^ y i3 t^kg^ ^ be equal to Tdust and 
chemical environment with which to probe their structure, excitation, radiative transfer and fluxes of the vari- 
properties, and evolution. Beyond their own intrinsic in- jnolecular lines are calculated for an assumed abun- 
terest, these regions may provide a hnk to the so-called ^^^^^ p^^gj^ ^^^^ ^ n(X)(r)/n(H2)(r). This trial abun- 
hot cores observed toward many high^mass star-forming ^^^^^ p^^^^^ ^^^^ ^^^.^^ ^^^^ ^^^^ agreement with 
regions (e.g., Walmsley & SchilkelM- observations is obtained. In practice, only two types of 
Rapid advances m both observational and modeling abundance profiles are considered: a constant abundance 
capabilities allow much more quantitative studies of the throughout the envelope or a 'jump' profile in which the 
chemistry in YSO envelopes than was possible even a few abundance is increased by a large factor in the inner warm 
years ago. Several different steps can be distinguished (see ^^^^^^^ jue to ice evaporation. Such models have success- 
Fig. Thanks to the advent of large-format bolometer f^^y ^^^^ applied to both high- (e.g., van der Ta k et al. 
arrays, most studies nowadays start with an analysis of the ^nd low-mass YSOs (e.g., CeccareUi et al. IMbl 

~ Z ~ ~ ~~ o ^ Schoier et al. 120021 . and work best for molecules for which 
bend offprint requests to: b. Doty ' 
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Fig. 1. An overview of methods used for constraining the physical and chemical structure of YSO envelopes. In general, 
the steps include the determination of the physical structure through dust modeling, calculation of gas temperature, 
and adoption of a chemical model. This combination produces observables (column densities, line fluxes, etc.) to 
compare with observations. The source parameters are determined by adjusting them until a best fit is obtained. Two 
important points of divergence include the use of self-consistent vs. approximate radiative transfer, and the use of 
a time-dependent chemical network ('Full chemical model') vs. simple trial abundances ('Empirical model') (figure 
adapted from van Dishoeck & van der Tak 1200(11 and van Dishoeck l2()()H|l . 



a large set of lines originating from levels with a range 
in energy has been observed. This method provides abun- 
dances for comparison with chemical models, but does not 
test the chemical networks directly. 

The second, ah initio or 'forward' approach is the 'full 
chemical model', in which only the density structure de- 
rived from the dust is adopted as a starting point. Given 
initial abundances and a detailed chemical network, the 
abundances of various molecules can be solved as func- 
tions of position and time and the gas temperature can be 
calculated explicitly by solving the full thermal balance of 
the gas. The physical structure n(r) and Tir) can either 
be taken to be constant with time or to vary according 
to some (dynamical) prescription. Such time- and space- 
dependent chemical models have been applied to low-mass 
YSOs by Ceccarelli, Hollenbach, & Tielens l(TU^ here- 
after CHT96) and Rodgers & Charnley and to 
specific high-mass sources by Millar, MacDonald & Gibb 



ifTMTI G34.3-H0.15) and Doty et al. PKI^ AFGL 2591). 
The main quantities to be determined are the best-fit time 
(or, in a dynamical model, mass-accretion rate) and other 
parameters which enter the chemical models, such as the 
cosmic ray ionization rate. The chemical models them- 
selves can be tested by comparing the abundance profiles 
ifit) at the best-fit time with those derived through the 
empirical method. In this way, they also provide a guide- 
line for more complicated abundance profiles to adopt in 
the empirical method. 

In this paper, we describe a detailed 'full chemical 
model' of the best studied low-mass YSO, IRAS 16293- 
2422. A novel feature is the addition of a full Monte Carlo 
radiative transfer calculation of the resulting line fluxes 
for direct comparison with observations (see bottom-right 
part of Fig. 0. Such models provide the most complete 
test of our understanding of the physical and chemical 
structure of YSO envelopes. 
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IRAS 16293-2422 is a nearby (- 160 pc, Whittet lWT^ 
low mass, low luminosity (27 L0) protostellar object lo- 
cated within the p Ophiuchus molcular cloud. It has an 
exceptionally rich and well-studied spectrum (e.g., Blake 
et al. 119941 van Dishoeck et al. 119951 Cecarelli et al. 
PO'OOa lOOOb, Schoier et al. lMT^ . and is therefore consid- 
ered the prototypical low-mass source for chemical stud- 
ies, much like Orion is for high-mass objects. Ceccarelli et 
al. (2000a, b) used the physical structure based on CHT96 
combined with a restricted chemical network to analyze 
data of H2CO, H2O, and SiO in a 'full chemical model', 
and found strong evidence for large abundance enhance- 
ments of these species in the innermost part (<150 AU) 
of the envelope. The evaporated species may subsequently 
drive a complex 'hot core' chemistry leading to the even 
more complex organic molecules which have recently been 
detected in IRAS16293-2422 (Cazaux et al. lMI^ . 

In a later analysis, Schoier et al. H2002|l combined 
dust/SED modeling of the physical structure, multiple 
line observations covering a range of excitation conditions, 
and a detailed radiative transfer analysis in an 'empir- 
ical model' to infer the structure of IRAS 16293-2422. 
This work supported the general conclusion of a 'hot core' 
where the abundances of key molecules are enhanced by 
several orders of magnitude due to evaporation of ices. 
The model employed only uniform and 'jump' abundances, 
however, which may not be representative of the detailed 
time- and space-dependent chemistry. 

Recently, Doty et al. 1)2002(1 described such a time- and 
space-dependent physical/chemical model for static YSO 
envelopes including the hot core chemistry. By combining 
the model results with observations of many species of 
one particular high-mass YSO, AFGL 2591, it has been 
shown that it may be possible to not only confirm the gross 
source structure, but also constrain source properties such 
as age, ionization rate, and role of grains in determining 
the chemical structure (see also Boonman et al. I2003|l . 
Here the 'full chemical model' of Fig. ^ was adopted, but 
the self-consistent line radiative transfer was performed 
for only a subset of the species. 

In this paper, we report on the application of the phys- 
ical/chemical model of Doty et al. I)2002|l to the low-mass 
YSO IRAS 16293-2422. These results are combined with 
a self-consistent radiative transfer model, and applied to 
the full multi-species, multi-transition dataset of Schoier 
et al. I(2002jl . By comparison with the case of AFGL 2591, 
we can also directly determine the differences in derived 
model parameters for a low- and a high-mass YSO (van 
Dishoeck 2003). The models and observations are briefly 
described in Section 2. The observations are then used 
with the models to constrain the source properties in 
Section 3. Finally, we summarize the results and conclude 
in Section 4. 



Table 1. IRAS 16293-2422 physical structure from dust 
modeling by Schoier et al. |2P02J) 



Parameter Value 



Distance, d (pc) 


160 




Luminosity, L (Lq) 


27 




Optical depth at 100 ^m, rioo 


4.5 




Density power law index, p 


1.7 




Inner envelope radius, (cm) 


4.8 X 


10^4 


Outer envelope radius, re (cm) 


1.2 X 


10^^ 


H2 density at 1000 AU, no (cm^^) 


6.7 X 


10^ 



2. Existing observations and models 

2.1. Observations 

IRAS 16293-2422 has been welLobserved both in the 
continuum and in various submillimeter molecular lines. 
While no new observations are presented in this paper, 
we briefly note and discuss the observational data as they 
provide the constraints placed on the model. 

The SED of IRAS 16293-2422 in the range 60 - 2900 
fim is presented by Schoier et al. (2002). High angular 
resolution JCMT data allowed the determination of radial 
brightness distributions at 450 and 850 fj,m on scales of 9" 
and 15" (1400 and 2400 AU) respectively. 

The molecular line data utilized here are taken primar- 
ily from the large surveys of IRAS 16293-2422 by Blake et 
al. l(1994|l and van Dishoeck et al. ljl995|l . Additional com- 
plementary data are taken from the JCMT public archive 
(Schoier et al..2002). These data are supplemented by the 
H2CO hues of Loinard et al. H2000|l . The data set - a to- 
tal of 76 transitions for 18 species considered here - has 
the advantage that it samples the full radial range of the 
envelope, providing probes over a wide range of physi- 
cal, thermal, and chemical conditions. Only information 
on the lowest transitions of the molecules, which occur at 
millimeter wavelengths and probe the very coldest outer 
parts and surrounding cloud, is lacking. In general, the 
calibration uncertainty of each individual line is ~30%. 

2.2. Model 

Here a brief synopsis of the physical, thermal, chemical, 
and radiative transfer models are provided. For more de- 
tailed information, see Doty et al. (2002), Schoier et al. 
(EOni, Doty & Neufeld ijTMTjl . and references there. 

2.2.1. Physical and thermal structure 

We adopt a spherically symmetric static model of the ex- 
tended envelope of IRAS 16293-2422 surrounding the two 
protostellar sources in the center (Looney et al. l2000f) . The 
observational continuum data were combined and simul- 
taneously modeled by Schoier et al. (2002) with the pub- 
licly available radiative transfer code DUSTY (Ivezic et 
al. 11999(1 . This analysis allows the source structure prop- 
erties (e.g., envelope size, density power law, continuum 
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Log [ Radial Position / cm ] 



Fig. 2. Physical and thermal structure of IRAS 16293- 
2422. The density and temperature are from the model of 
Schoier et al. H2(JU2fl . The gas temperature is assumed to 
follow the dust temperature. 



optical depth) to be determined to within approximately 
±20% (see also Doty & Palotti lMTl . The adopted source 
properties are presented in Tabled and the density and 
temperature structure are reproduced in Fig. El 

The successful line modeling of Ceccarelli et al. I)2()()()a|l 
and Schoier et al. (2002) in spherical symmetry down to 
~ 30 AU suggests that the assumption of spherical sym- 
metry may be largely justifiable. Most chemical data were 
obtained with a 15" beam, which yields a linear size of 
2400 AU at the assumed distance of 160 pc, larger than 
the 800 AU separation of the central protostars. Schoier 
et al. H2003|l recently used interferometric obesrvations to 
determine the structure below 1000 AU, finding that the 
binary has cleared out most of the material in the inner 
part of the envelope, and that there exist two unresolved 
central sources with best-fit disk sizes of ~ 250 AU in 
diameter. However, they find that the while the detailed 
results for the inner envelope leave the inner (r < 400 AU) 
structure somewhat uncertain, their results have little ef- 
fect on the extended envelope (r > 400 AU). 

A majority of the modeled envelope exists at relatively 
high densities and optical depths, leading to a strong ther- 
mal coupling between the gas and dust. As a result, the 
gas temperature is assumed to follow the dust temperature 
(CeccareUi et al. lTO^ Doty & Neufeld lTWTI CeccarelU et 
al. l2000a|l . Test calculations have shown that this assump- 
tion is sufficient for both the chemistry (Doty et al. l2002|l 
and radiative transfer through molecular lines (Boonman 
et al. lMIHji . 



2.2.2. Chemistry 

The chemical model is based upon the UMIST gas-phase 
chemical reaction network (Millar, Farquhar, & Willacy 
119971 hereafter MFW), including reactions to model the 
hot-core chemistry. Pseudo time-dependent models of the 
chemical evolution over 30 radial grid points were con- 
structed, providing a time- and space-dependent chem- 
ical evolution. The local parameters (hydrogen density, 
temperature, and optical depth) at each radial point are 
taken from the physical and thermal structure calculations 
above. For the majority of species, the initial abundances 
of the high-mass source AFGL 2591 (Doty et al. l2002|l . are 
utilized as shown in Table |2] The chemistry of deuterated 
molecules is not considered. 

The effects of freeze out onto and desorption from dust 
grains are approximated. Instead of explicit freeze out or 
desorption with time, the desorption is taken to be instan- 
taneous at 100 K, where expected grain mantle species are 
injected into the gas fCharnlev [ig?7|l . in keeping with the 
timescales observed in the laboratory (Fraser et al. l2001|l . 
Below 100 K, we deplete the gas phase abundance of many 
species expected to be in ices (e.g., H2O). The only excep- 
tions are CO, which we take to desorb at Tco, and II2CO 
and CH3OH, which we take to desorb or 'jump' at Tx- 

The effects of photodissociation from the interstellar 
radiation field at the outer boundary is included, but is 
generally small due to the high optical depth, and the 
coarseness of the spatial grid considered. Finally, we have 
considered the effect of metal depletion by significantly 
reducing the initial Fe abundance. We find that metal de- 
pletion makes only a small difference, worsening the fits 
by only a few percent on average. 

2.2.3. Radiative transfer 

The molecular line radiative transfer is accomplished 
through a non-LTE, Monte-Carlo model described in 
Schoier l|2()()()|l . This code has been benchmarked to high 
accuracy against a wide range of other molecular line ra- 
diative transfer models (van Zadelhoff et al. 12002(1 . In this 
model, the spatial molecular abundances x{r,t) are com- 
bined with the adopted physical structure to compute the 
excitation and resulting line profiles for all transitions up 
to ~ 500 K in the ground vibrational state of the observed 
molecules. Chemical evolution times from 3 x 10^ years to 
3 X 10^ years are considered, with one dex spacing. 



3. Results 

In this section the results of our physical/thermal/ chem- 
ical modeling of IRAS 16293-2422 are presented, and the 
comparison of line strengths predicted from this model to 
those observed. As a metric of the comparison, we adopt 
the mean percentage magnitude difference between the 
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Table 2. Initial abundances at t—0 relative to H2 
AFGL 2591 



for 



Species 


Initial Abundance 


Ref. 


Initial Abundance (T > lOOK) 






CO 


3.7(-4) 


a 


CO2 


3.0(-5) 


d 


H2O 


1.5(-4) 


d 


H2S 


1.6(-6) 


h 


N2 


7.0(-5) 


e 


CH4 


1.0(-7) 


e 


C2II4 


8.0(-8) 


e 


C2II6 


1.0(-8) 


e 





0.0(0) 


e 


H2C0 


1.2(-7) 


e 


CH30H 


1.0(-6) 


e 


s 


0.0(0) 


e 


Fe 


2.0(-8) 


e 


Initial Abundance (T < lOOK) 






CO 


3.7(-4) 


a 


CO2 


0.0(0) 


f 


H2O 


0.0(0) 


f 


H2S 


0.0(0) 


f 


N2 


7.0(-5) 


e 


CH4 


1.0(-7) 


e 


C2II4 


8.0(-8) 


e 


C2II6 


1.0(-8) 


e 





8.0(-5) 


g 


H2CO 


0.0(0) 


f 


CH3OH 


0.0(0) 


f 


s 


6.0(-9) 


h 


Fe 


2.0(-8) 


e 



a{b) means a x lO'' 
All abundances are gas-phase, and relative to H2 
van der Tak et al. T999 ^ van der Tak et al. l^nUTJl " van der 
Tak & van Dishoeck .20Q0. Boonman et al. l^ITOl " Charnley 
119971 ^ assumed frozen-out or absent in cold gas-phase, ^ taken 
to be ~ consistent with Meyer, Jura, & Cardelli 1998 ^ Doty 

et al. i^nn^ 



predicted and observed line strengths, given by 



A = 



mod.i 



obs.i 



obs.i 



(1) 



This form has the advantage that it measures the size 
of the difference, without allowing equally balanced high 
and low values to cancel out. Furthermore, the summation 
over lines implicitly provides a greater weight to molecules 
for which more data - and thus more constraints - exist. 
Measures using the sum over molecules instead of lines 
show qualitatively similar, though often accentuated, re- 
sults. We note that we adopt this measure instead of a 
analysis, as agreement to within a factor of 10 are con- 
sidered good for chemical modeling. This value is enough 
larger than the statistical uncertainty in the line obser- 
vations so as to invalidate the statistical meaning of 
measure. 

The parameters varied in the models are the cosmic 
ray ionization rate ^, the adopted initial abundances in the 
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Fig. 3. The ratio of the predicted and observed line 
strengths for molecules observed toward IRAS 16293- 
2422. The errorbars represent the range of values between 
3 X 10^ < i(yrs) < 3 x 10"^. This represents the best fit time 
range, and is consistent with previous values based upon 
estimates of the infall rate and central mass. Guidelines 
at factor of 3 (dotted) and 10 (dashed) ratios are given to 
indicate good and acceptable fits. 

inner and outer regions, and the desorption temperatures 
of selected species (CO, H2CO, CH3OH). Detailed radial 
profiles of selected species are presented in §3.7 and 3.8. 

3.1. General Results 

As discussed above, our base model is taken after the high- 
mass hot core + envelope of AFGL 2591 by Doty et al. 
f2002). A parameter space search, guided by results from 
previous studies, suggests a best fit to the observed data 
of IRAS 16293-2422 with only minor modifications to the 
AFGL 2591 initial chemical conditions. 

A comparison between the best fit model and observa- 
tions is shown in Fig.O Here the ratio of the predicted to 
observed line strengths for each molecule observed is plot- 
ted. We find a best-fit time of 3 x 10^ < t(yrs) < 3 x 10^, 
with these times forming the range in the figure. These 
times are consistent with the age inferred by Schoier et al. 
(|2002| from fitting a collapse model to the line profiles, 
and by using the constant infall rate of Ceccarelli et al. 
H2000a|l with their preferred central mass of 0.8 Mq. 

As can be seen in the figure, the majority of species (11 
of 18) are fit to within 50 % of the observations, thirteen 
arc fit to within a factor of three, and 15 are fit to within 
a factor of 10, a level usually considered acceptable agree- 
ment in chemical modeling (see e.g., Millar & Freeman 
11984^ Brown & Charnlev ITMH Terzieva & Herbst [MHI) ■ 
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Table 3. Differences between best fit model for IRAS 
16293-2422 and AFGL 2591 



Parameter 



IRAS 16293 AFGL 2591 



Initial Abundance (T > lOOK) 










CO 


1.0(- 


-4) 


3.7(- 


-4) 


HaS 


1.0(^ 


-8) 


1.6(- 


-6) 


H2CO 


8.0(- 


-8) 


1.2(- 


-7) 


CH3OH 


1.5(- 


-7) 


1.0(- 


-6) 


Initial Abundance (T < lOOK) 










CO (20 < T(K) < 100) 


1.0(- 


-4) 


3.7(- 


-4) 





1.0(- 


-4) 


8.0(- 


-5) 


H2CO (60 < T(K) < 100) 


8.0(- 


-8) 


0.0 




CH3OH (60 < r(K) < 100) 


1.5(- 


-7) 


0.0 




Tde.(K) [X{T < Tdes) ~ 0] 










CO 


20 




100 




H2C0 


60 




100 




CH30H 


60 




100 




CR ionziation rate (s^^) 


5.0(- 


-16) 


5.6(- 


-17) 



a{h) means a x 10*" 
All abundances are gas-phase, and relative to H2 



An interesting case is the ~ 20 — 30% deviation in ^^CO, 
and the small uncertainties on ^^CO, C^^O, and HCO+. 
It is possible that the ^^CO discrepancy could be due to 
a different ^^G/^^G ratio than taken here. It may also be 
due to deviation of the real structure from the continu- 
ous, spherical symmetry we have adopted. In any case, 
the deviation is no larger than the expected calibration 
uncertainty of ^ 20 — 30%. 

In the comparison there exist three species which de- 
viate by more than an order of magnitude. The outliers 
are: OCS, C2H, and HCS+, which have individual devia- 
tions of a factor 30, 30, and 100 respectively. While there 
is variation, our model tends toward producing too little 
emission. The difference between the model and observa- 
tions is A = 0.51 when these are omitted. Including them 
raises A to 0.87. This quantitatively confirms the agree- 
ment between the model and observations at the level of 
a factor of two for most species. 

The differences between the AFGL 2591 model and the 
best fit model here are summarized in Table 13 As can be 
seen, the dilferences are generally minor. They are: (1) a 
large increase in the cosmic-ray ionization rate by a factor 
of >10 over the 'standard' value of ^ 10"^^ s~^ and the 
AFGL 2591 value of - 5.6 x lO^i^ s^^; (2) depletion of CO 
at low temperatures 20 — 40 K); (3) depletion of H2CO 
and CH3OH at moderate temperature (< 60 K); and (4) 
variations in the initial abundances of a few other species. 
We discuss each of these differences below separately, as 
variations from the best fit model. 



0.7 



0.65 



0.6 



0.55 



0.5 - 



0.45 




-17 -16 -15 

Cosmic-Ray Ionization [(J'] (s~') 



-14 



Fig. 4. Dependence of quality of fit as measured by 
the mean difference between predicted and observed line 
strengths, as a function of the ionization rate, C,. Notice 
the best fit near ^ = 5 x 10~^^ s""'^. The two points near 
the minimum are indistinguishable at the level of uncer- 
tainty of the observations, and given the constraints of the 
model. In either case, the ionization rate is much larger 
than standard. 



3.2. Effects of cosmic ray ionization rate 



'16 



The ionization rate inferred in the modeling, ^ 5 x 10 
s^^, is much higher than the 'standard' cosmic ray ion- 
ization rate of ~ lO"^"^ s'^ (e.g., Roberts & Herbst 1 ^170^ 
Black k Dalgarno Unni O'Donnell & Watson il974p . To 
see the dependence of the model results on the ioniza- 
tion rate. Fig. ^ shows the mean difference A between 
the models and observations upon varying the ioniza- 
tion rate. The minimum deviation occurs in the range 
^ = 5 X 10~^^ — 10~^^ s~^. While the two values near the 
minimum are essentially indistinguishable, the ionization 
rate required for this fit is 50-100 times higher than the 
traditional cosmic-ray ionization rate used in dark cloud 
models (Lepp 1992). 

The best fit range for C is expected to be meaning- 
ful, due to the fact that the mean difference is minimized 
here. The variation in A is damped by the fact that it is 
an average across all species. As a result, a 25% variation 
in A can correspond to a factor of 2 change in 1/4 of the 
species, or a factor of ~ 8 change in ^ 4 — 5 transitions. 
This is seen for the case of CO in Figures [T] and |S1 where a 
smaller change in A over all species corresponds to vari- 
ation in physical parameters by ^ 100 x, and variations 
in line strengths by up to 8 x. As a result, we infer that 
a minimum in A and a variation of ^ 20% is sufficient to 
draw conclusions, which implies that the preferred value 
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of C = 5 X 10 ^® — 10 s ^ is meaningfully different from 
other values tested. 

The ionization rate C is assumed to be uniform 
throughout the source. Species that are most affected by 
the variation in the ionization rate are HCO+, HCN, SO, 
and H2CO and show improvements of up to 100%. While 
the 1100+ abundance should be directly related to the 
ionization rate throughout the envelope, Doty et al. : 2()()2} 
show that the remainder are predominantly active above 
100 K. This implies that the ionization rate may be po- 
sition dependent, with the most affected species in the 
interior. 

The origin of this enhanced ionization is of physical 
interest. Recent models and measurements infer up to 2 
orders of magnitude variation in the cosmic ray ionization 
rate (McCaU et al.'MJT Liszt MF: Doty et aLESSI van 
der Tak & van Dishoeck 2(1(1(1; . It is difficult to understand 
such an extreme variation from source to source, imply- 
ing that some other physical mechanism may produce or 
contribute to the ionization. This is especially true as the 
recent 'high' inferred values for ( are for diffuse clouds, 
while dense cloud models have historically required much 
smaller cosmic-ray ionization rates near 1 — 3 x 10~^^ s~^. 
In the case of IRAS 16293, we suggest the possibility that 
the enhanced ionization is due to X-rays produced by mag- 
netic activity associated with accretion onto the proto- 
stars. This could both produce the exceptionally high in- 
ferred ionization rate, and preferentially affect the warmer 
interior species. The effects of a central X-ray source will 
be presented in a forthcoming paper (Doty et al. , in prepa- 
ration). 

3.3. Effects of CO desorption temperature 

Studies of sohd CO on icy dust grains (Colhngs et al. l2UU3l 
Eraser et al. IMH Galloway & Herbst 1994; Sandford & 
Allamandola ll993al Nair & Adamson ll97();i suggest that 
the bulk of the CO evaporates in a step-wise fashion be- 
tween 20 and 70 K - depending upon whether it is trapped 
inside or lies on top of the ice - much lower than the tem- 
perature at which water ice desorbs from grains (Fraser 
et al. I2()()l|l . This is consistent with our best fit model. 
In order to test this, we have varied the CO desorption 
temperature from our baseline best-fit model. The results 
are shown in Fig. |S1 

Clearly, the best fit requires a CO desorption temper- 
ature near ^ 20 K and < 60 K. A lower temperature 
both yields a worse fit to the observational data by over- 
producing the ^'^CO and C^^O emission by 25% and 68% 
respectively at Tqq = 10 K, and is inconsistent with lab- 
oratory results. Much higher temperatures yield signifi- 
cantly worse fits to the data, underproducing both ^'^CO 
and C^^O hue fluxes by a factor of 50 by Too = 100 K. 
The species most affected (aside from CO itself) are the 
cyanogens CN, HCN and HNC, and the CO ion-molecule 
byproducts HCO+, CS, and II2CO, all of which show vari- 
ations between 40-300% . 
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Fig. 5. Dependence of quality of fit as measured by 
the mean difference between predicted and observed line 
strengths for the full chemistry / observational set, as a 
function of the CO desorption temperature (Too)- Here, 
j:(CO) = 10^4. Notice the best fit for Tco ~ 20 K. 



Physically, a low desorption temperature near 20 K 
would be an indication that a significant fraction of the 
CO is not intermixed with the H2O in the grain mantle. 
A number of suggestions for differentiation in the ice have 
been made, including differentiation in the gas prior to 
adsorption, differentiation in the ice due to chemical and 
physical processing, and differentiated freeze-out during 
the cooling time behind a shock which has liberated the 
grain mantles (e.g., Schutte [19971 ; Bergin et al. 11999(1 . 
While it is difficult to comment on the first two scenar- 
ios, it is doubtful that shock processing is the main cause 
for IRAS 16293-2422. In particular, the products of shock 
chemistry do not dominate the bulk of the envelope, and 
the small lincwidths observed for many species further 
suggest that a large fraction of the volume of the gas is 
not shocked. Note that this analysis does not exclude that 
some fraction of the CO also evaporates at higher tem- 
peratures. In fact, J0rgensen et al. f2002j) conclude from 
their analysis of the CO 3^2 and 2^1 isotopic lines 
in a sample of class objects that some CO must still be 
frozen out at temperatures above 25 K. 

3.4. Effects of CO abundance 

Previous authors (Schoier et al. 120021 Ceccarelli et al. 
I2000a|l have inferred (constant) CO abundances in the 
range of 10~^ - 10""'. As a result the inner (T > Tdcs,co) 
CO abundance is varied, keeping all other parameters the 
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Fig. 6. Dependence of quality of fit as measured by 
the mean difference between predicted and observed fine 
strengths as a function of the CO fractional abundance for 
T>Tco= 20 K. Notice the best fit for a;(CO) - IQ-". 

same as in our best fit model. The results are presented 
in Fig. ini 

As can be seen, CO abundances of 1 — 4 x 10~* are 
preferred. Lower abundances produce too little C^^O emis- 
sion to be consistent with observations. The best fit ~ 
based upon CO data alone - is a;(CO) — 10"'* where mod- 
els match C^^O observations to within 2%. By a;(CO) = 
5 X 10^^ and 4 x 10^*, the discrepancy for the C^^O 
lines reaches 37% and 95% respectively. The effect on A is 
smaller for two reasons: first other species are included in 
the mean difference, and second the comparison for HCO+ 
and CN are both improved as x{CO) increases. 

The HCO^ and CN abundances both rely upon CO via 
ion-molecule reactions. Since this relation is much more 
indirect, the roles of HCO+ and CN abundances in fixing 
the CO abundance are discounted, and a:;(CO) ^ 10^* is 
preferred. The constant abundance of ^ 4 x 10^^ inferred 
by Schoier et al. PIU^ in the empirical model can now 
be understood as a weighted average of a very low CO 
abundance at T < 20 K and a higher abundance of ~ lO^'* 
at T > 20 K. 

We have also considerd the combined effects of CO 
desorption temperature, cold (T < Tqo) CO abundance, 
and warm (T > Tqq) CO abundance. The results are 
presented in Fig. \7\ where we plot the mean percentage 
difference between model and observed line strengths for 
C^^O (top panel), and all CO isotopomers (bottom panel). 
Warm CO abundances are denoted by the different line 
types, corresponding to a::(CO, warm) = 10^"*, 5 x 10^^, 2 x 
10~^ respectively. As can be seen, the CO data strongly 
prefer a;(CO,warm) ~ 10"'', consistent with our results 
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Fig. 7. Contours of the percentage difference between 
model and observed CO line strengths, for various values 
of CO desorption temperature (Tco), abundance of cold 
CO (T < Tco), and abundance of warm CO (T > Tco) 
denoted by the different line types. The solid, dotted, and 
dashed lines correspond to x{CO, warm) = 10"'*, 5 x lO"'*, 
and 2 x 10~^ respectively. The top panel gives the re- 
sults for C*^0, and the bottom panel the average results 
for i^CO, 0*^0, and 0*^0. Notice that the 0*^0 results 
strongly favor a;(CO,warm) — 10^*, as do the results av- 
eraged over the CO isotopomers. 

above. This is confirmed by the results in Fig. |S1 where 
we plot contours of the mean difference. A, over all ob- 
served lines. Again, models with a;(CO, warm) ~ 10~^ are 
preferred. 

Perhaps even more interestingly, the results for both 
the CO and general chemical network allow us to simul- 
taneously constrain the desorption temperature and cold 
CO abundance. Taking uncertainties in the observational 
data of ±30% suggests 18 < Tco(K) < 23 based upon the 
25% contour level in Fig. [7| Even if the uncertainties are 
larger, the results in Figs. [7| and |H1 provide outer limits of 
16 < Tco(K) < 30. There may be a potential region of de- 
generacy in Figs.[7|&|Hlas the cold CO abundance appears 
to increase with increasing Tco- It should be noted that 
the cold CO abundance is still ^ 10~^ in this result 
that may be explained by significant evaporation near Tco 
and some partial/gradual evaporation of CO presumably 
in an H2O matrix at higher temperatures. These results 
are consistent with both the cut along Tco and the labo- 
ratory results discussed above. 

Finally, there does seem to be evidence of depletion 
in the cold gas for T < Tco- The comparison for both 
the CO and overall set of observed species suggest a rela- 
tively firm upper limit of a; (CO, cold) < 10~^. The regions 
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CO Desorption Temperature (K) 

Fig. 8. Contours of quality of fit as measured by the mean 
difference between predicted and observed hne strengths 
over all species as a function of the CO desorption tem- 
perature (Tco), abundance of cold CO (T < Tco), and 
abundance of warm CO {T > Tco) denoted by the differ- 
ent line types. The solid, dotted, and dashed lines corre- 
spond to 2;(C0,warm) = 10"", 5 x 10"'', and 2 x 10"^ 
respectively. Notice the agreement with Fig. in con- 
straining 16 < Tco < 30 K, a;(CO,cold) < 10-^ and 
x(CO,warm) - 10"'*. 



of best fit appear to encompass values of 3-30 times less 
(3 X lO"'^ to 3 X 10^6). The upper value of 10"^ signifies 
a depletion of 90%, while the lower values correspond to 
97% and 99% depletion respectively. While the exact level 
of depletion is uncertain, these results do confirm a sig- 
nificant sink of gas-phase CO - presumably as ices onto 
dust grains in the cool exterior. Such high levels of CO de- 
pletion are consistent with those found in cold pre-stellar 
cores (e.g., Bacmann et al. l2003|l and the large abundances 
of deuterated molecules detected in the outer envelope of 
IRAS 16293-2422 (van Dishoeck et al. lTM^ Loinard et al. 
IMTHI Parise et al. IMI^ 

3.5. Effects ofH2C0 and CH^OH depletion 
temperature 

We also consider the effects of modifying the H2CO 
and CH3OH desorption temperatures, Tx- In the case of 
II2CO which may be produced at some level in the gas 
phase (Doty et al. I2002| l , this corresponds to the temper- 
ature at which significant production occurs. The results 
are shown in Fig. where the mean difference is plotted 
as a function of Tx- 
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Fig. 9. Dependence of quality of fit as measured by 
the mean difference between predicted and observed line 
strengths as a function of the temperature of II2CO and 
CH3OH desorption [Tx). Notice the best fit for Tx > 30 
K. 

To within the uncertainties of the observations and so- 
phistication of the models, the results for Tx = 30 — 100 
K are indistinguishable. Consequently, the results in Fig. 
El suggest Tx > 30 K. In the case of CH3OH, which is 
almost certainly formed on the grain surface (e.g., Tielens 

6 HagenHlHl Blake et al. DESL Doty et al. lM)^ . this 
is probably an indication of desorption. Desorption in 
this temperature range is consistent with current chem- 
ical understanding of solid ices. Sandford & AUamandola 
Ill993bll measure a pure CH3OH desorption temperature 
(in space) of 70 — 80 K. This is in keeping with the fact that 
CH3OH should have a lower desorption temperature than 
water due to its weaker hydrogen bonding. For compar- 
ison, a Clausius-Clapeyron calculation which reproduces 
the water evaporation temperature well suggests an evap- 
oration temperature for CH3OH of 87 K (Alsindi et 
al. 12003(1 . Likewise, these results are consistent with the 
empirical modeling of Schoier et al. f 2002,1 . who found 
Tx(CH30H) > 50 K. As a result, it is encouraging that 
laboratory work, theoretical calculations, empirical mod- 
els, and the results of this work are all in agreement with 
a CH3OH desorption temperature of 60 < r(K) < 100. 

In the case of H2CO, the meaning of this tempera- 
ture is less clear. The deviation between the observed and 
predicted II2CO line strengths caused by the overproduc- 
tion of II2CO in the model do not significantly change 
as Tx is increased - from 49% at T^cs = 10 K to 44% 
at Tdcs = 100 K. While Doty et al. ijMT^ suggested it 
was possible that gas-phase reactions could play a role in 
the H2CO chemistry, they did not identify any that would 
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cause a significant 'jump' in abundances in this temper- 
ature range. While it is difficult to directly constrain the 
H2CO abundance, it is clear that the modeling here is in- 
sensitive to the amount of cold H2CO, and that while no 
jump is required, a jump due to desorption is not ruled 
out - consistent with Tx > 40K as found by Schoier et al. 

mm . 

3.6. Outliers 

As can be seen in Fig. |21 the three species OCS, C2H, 
and HCS"*" are outliers in our models. These species yield 
line strengths which diverge from the observations by fac- 
tors of '-^ 5, 100, 30, and 100 respectively. Some deviation 
due to radiative transfer, geometrical, and line of sight ef- 
fects is to be expected. However, the differences for these 
three species are discrepant from the other species consid- 
ered, implying that while our adopted physical/chemical 
structure is reasonable, our knowledge of the chemistry is 
lacking. 

In the case of OCS, the chemistry is uncertain, and 
many of the reaction rates are estimates without signifi- 
cant laboratory study. As a result, it is not suprising that 
a discrepancy exists. For C2H, UV photodissociation in 
the outer region - while included - is not significant as 
most of the C2H is produced above lOOK. The dominant 
production is via recombination of C2H5", and reaction of 
CsH^ with O. Destruction is mainly through the neutral- 
neutral reaction with O. While a lower oxygen abundance 
can increase the C2H abundance, it also decreases the 
abundances of the other important oxygen-bearing species 
such as SO, SO2, and CH3OH to below the observations. 
On the other hand, it is interesting to note that the de- 
struction reaction with O is assumed to be temperature 
independent (MFW). The existence of a reaction barrier 
or a temperature dependent rate of collisions would both 
tend to decrease the destruction, which would have the 
effect of raising the C2H abundance closer to the observed 
levels. 

The third outlier, HCS"*" follows an ionization balance 
with CS via dissociative recombination, and reactions with 
HCO+, H+, and H3O+. However, raising ( to 5x lO^^^ s^^ 
only changes the HCS+ discrepancy by ^ 3%. Such high 
values of the cosmic ray rate raise the most sensitive ion - 
HCO"*" - to levels far above that observed. While the ma- 
jor destruction mechanism of dissociative recombination 
has a somewhat strong temperature dependence, 

this rate has been measured in the laboratory, and is con- 
sidered to be accurate to within 25% by MFW. This leaves 
the production reactions of (HCO"'", Hj^, and HsO"*") -I- CS 
as potential sources of uncertainty. Each of these rates are 
estimated. As such, it would be useful to measure them in 
the laboratory to confirm the rates adopted by MFW. 

Finally, while HNC only differs by a factor of 5 and 
is thus not a major outlier, the production of HNC is 
still not fully understood (e.g. Rodgers & Charnlev 120011 
Charnley, Rodgers & Ehrenfreund 120011 Liszt & Lucas 



12001(1 . In our model, HNC is produced primarily through 
the dissociative recombination of HCNH+ and H2NC+. 
In the recombination of HCNH+, where the branching 
fractions are taken to be 50% for CN, 25% for HON, 
and 25% for HNC, the HNCH+ abundance is determined 
mostly by an ionization equilibrium in which the primary 
production paths are proton transfer between HNC and 
(HCO^, H3O+, H;^). These are also the dominant de- 
struction paths for HNC. On the other hand, H2NC'^ 
is formed by C+-f NH3 H2NC'^, and dissociatively 
recombines to form HNC and CN in a 10:1 ratio. We 
note that while too much HNC is produced in our model, 
the CN abundance is somewhat low. This combination 
suggests that the adopted branching fractions for the 
dissociative recombination should perhaps be reinvesti- 
gated. Adopting the results of Talbi & Herbst H1998|l for 
C+ + NH3 HCNH+ has little effect, suggesting fur- 
ther concentration on H2NC+ Alternatively, some HNC 
destruction route may be missing in the networks. 

3.7. Comparison with empirical model 

In principle, the empirical modeling approach should - 
with sufficient parameter variation - mimic the results of 
the full chemical modeling. In practice, it is difficult to 
parametrically vary all species in a sufficiently meaning- 
ful yet complete manner. As discussed in §1, the approach 
adopted by many authors is to treat abundances as ei- 
ther constant, or as piecewise constant with 'jumps' at 
appropriate temperatures. This was the approach taken 
in Schoier et al (2002 ). It is interesting to compare the 
inferred abundance distributions from the empirical mod- 
eling from the more detailed results of the full chemical 
modeling. 

Quantitatively, the Schoier et al. <|2002|l empirical 
abundances reproduce the observations with A = 0.24, 
about half of our best fit model, A = 0.51. This is, how- 
ever, not a suprise. In the Schoier et al. H2002|l model, 
there are many more free parameters as the abundances 
for each of the observed species are varied both above and 
below the assumed desorption temperature. Furthermore, 
these variations are done without respect to constraints 
on the chemical network or evolution time. On the other 
hand, we directly specify the initial abundances for only 
three of the observed species (CO, H2CO, and CH3OH), 
and are constrained by the chemical network and its evo- 
lution. Furthermore, as discussed previously, the majority 
of the inputs to the chemical network are taken directly 
from agreement reached on a high-mass hot-core source, 
AFGL 2591. Taken together, these results are strongly 
encouraging as the chemical network comparison gives a 
good fit for significantly fewer direct parameter variations, 
confirms the proposed age of the source from physical evo- 
lution models, and directly tests the validity and extensi- 
bility of the chemical networks. 

As a more direct comparion, the radial abundance pro- 
files for the species discussed in this paper are presented in 
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Fig. 10. Radial abundance profiles for the species considered in the text for t = 3 x lO'' and 3 x 10'' years. The results 
of the empirical model of Schoier et al. H2UU2|I are given by the lines with arrows, for comparison. 



Fig- El where the two limiting times of 3 x 10'^ years and 
3 X 10^ years are plotted. In general, the abundances in- 
ferred from the empirical modeling are grossly consistent 
with those from the more detailed full chemical modeling. 
For most species with constant abundances, the inferred 
abundances in mostly the outer envelope are equivalent 
to within a factor of 3-10 at our best-fit time (here taken 
to be that intermediate between the two limits). In some 
cases, e.g. CN, CS, and CH2CO, the chemical model abun- 
dance oscillates with position in the cloud and the empir- 
ical value is simply a rough average of these complicated 



profiles. The most significant discrepancies are for OCS, 
HCS''",C2H, and HNC, which are discussed above. 

Even more interesting is the comparison of the jump 
models. In general, those species which show significant 
spatial variation in the full chemical model are represented 
as 'jumps' in the Schoier et al. (2002!| empirical model. Of 
these, the empirical and full chemical models generally 
agree to within a factor of 3 or so. As discussed in §3.4, 
the inferred CO abundance can be understood in such a 
'jump' model. The significant discrepancy is H2CO, which 
the full chemical model predicts to have only a small jump, 
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while the empirical model infers a cold, outer abundance 
some 15 times lower. 



3.8. Predictions for future observations 

While the model presented is able to simultaneously match 
the SED and many of the molecular line observations, a 
significant test will be the predictions it makes that can be 
studied by future facilities. In particular, with CARMA, 
SIRTF, and SOFIA data to be available in the next few 
years, and the upcoming leaps in resolution and sensitivity 
from Herschel and ALMA, it will be possible to probe 
many of the transitions and much of the spatial structure 
of IRAS 16293. 

To aid such future observations, the radial distribu- 
tions of a number of interesting species are shown in Fig. 
111! Also, column densities predicted by the model toward 
IRAS 16293-2422 for t = 3 x 10^ years and 3 x 10^ years 
are given in Table 01 This includes, in particular, radial 
column densities given by iVx.rad = J nx{r)dr, where 
nx{r) is the density of species AT as a function of r. We 
also include the column density averaged over a 15 arc- 
second beam toward IRAS 16293, given by A'^x.boam — 
J j nx{z,b)dzG{b)2TTbdh/ J G{b)2Tihdb, where b is the im- 
pact parameter, and G{b) is the beam response function. 
The column densities are sorted from highest to lowest in 
the 15 arcsecond beam at t = 3 x 10^ years, and continue 
down to ^ IQ^^ cm^^. 

Finally, it is intruiging to speculate on the effects and 
potential observability of a collapsing cloud. The evolution 
of a gas parcel in a protostellar envelope passing from a 
low temperature / density to a high temperature / density 
through infall has been studied by CIIT96, and Rodgers & 
Charnley l)2()()3|l among others. In the outer regions, cool 
ion-molecule chemistry dominates. As the gas heats while 
infalling, each adsorbed species passes through a sublima- 
tion front. In the interior warm neutral chemistry can take 
place, for example the conversion of gas-phase oxygen to 
H2O (e.g., CHT96; Charnlev ITM7|l . and the production 
of significant complex cyanogens and hydrocarbons (e.g., 
Rodgers & Charnley "^007, Doty et al. In free-fall 

collapse of low-mass YSOs, the dynamical timescale in the 
warm interior is less than the chemical timescale, leading 
to inner core abundances that mirror those in the cool 
exterior. The observations of a wide variety of complex 
daughter species observed in these warm interiors (e.g. 
Cazaux et al. I2()()3|l can only be understood if the col- 
lapse/infall is "slowed" to timescales over at least 10** 
years (Rodgers & Charnley I2003|l - a scenario more in- 
line with our adopted static case. Since the chemistry en- 
codes the temperature/density temporal evolution of the 
gas, small-scale spatial variations in abundances should, in 
principle, be able to distinguish between various dynami- 
cal scenarios such as static, Shu (,1977,) collapse, Larson- 
Penston (LarsonHHSi Penston HM^ infall, etc. This will, 
however, require a next generation of modeling that in- 
cludes detailed dynamics, thermal balance, chemistry, & 
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Fig. 11. Radial abundance profiles for the some of the more abundant species that may be targets for future observa- 
tions. The data are reproduced for i = 3 x 10'^ and 3 x 10** years. 



radiative transfer, as well as high spatial and spectral res- 
olution observations with instruments such as Herschel, 
ALMA, CARMA, and SOFIA to probe multiple lines at 
^ 100 AU resolution. 

4. Conclusions 

We have constructed detailed thermal and gas-phase 
chemical models for IRAS 16293-2422 based upon the 
physical model of Schoier et al. H2()02|l . These models 
were used to probe the validity of the proposed physi- 



cal structure, as well as study the chemical evolution of 
the source, and to test the application of our combined 
'hot-core'/envelope chemistry model of AFGL 2591 to a 
low-mass 'hot-core'-like source. In particular, we find that: 

1. The combined application of a physical, thermal, and 
chemical model with detailed radiative transfer is a 
powerful tool in constraining the structure and evolu- 
tion of depth-dependent sources. 

2. The time and position dependent model of Doty et al. 
(.2002,) can be meaningfully applied to IRAS 16293- 



S. D. Doty et al.: Chemistry in IRAS 16293-2422 



14 

2422, yielding results qualitatively similar to the mas- 
sive YSO AFGL 2591. 

3. The best fit for IRAS 16293-2422 occurs for times in 
the range 3 x lO'^ < t(yrs) < 3 x 10'*, consistent with 
existing infall models (Sect. 3.1). 

4. The best-fit ionization rate in IRAS 16293-2422 is high, 
~ 5 X 10~*^ — 10~*^ s~*, compared with previous re- 
sults for dense clouds. We propose that this may be 
due to X-ray emission from the central sources (Sect. 
3.2). 

5. Our best fit suggests that important CO desorption 
occurs at low temperatures, ~ 20 K, and constrained 
to the range 15 < T'co(K) < 40. Some sohd CO may 
remain at higher temperatures for these timescales. 
These results are in agreement with recent laboratory 
data of CO on - but not mixed with - a water ice 
(Sect. 3.3). 

6. We can also constrain the warm (T > Tco) and cold 
(T < Too) CO abundances to a:(CO,warm) - lO"-* 
and a:(CO,cold) < 10^^ and most probably < 10"^. 
These results are reflected in both the CO lines and 
in the results from the greater chemical network, and 
suggest significant (> 90%) depletions at low temper- 
atures (Sect. 3.4). 

7. CH3OH appears to desorb at temperatures 60 < 
T(K) < 100, consistent with laboratory expectations. 
On the other hand, the comparison between the H2CO 
predictions and observations are insensitive to the 
amount of cold H2CO present (Sect. 3.5). 

8. The chemistry of HNC, C2H, and HCS+ may not be 
fully understood. In particular, it would be useful to 
measure the branching fraction of dissociative recombi- 
nation of H2NC"'", the temperature dependence of the 
reaction 0H-C2H, and the ion-molecule reaction rates 
(HCO+, H+, H3O+) + CS (Sect. 3.6). 
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